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ABSTRACT 

Godunov's method and several other methods for computing solutions to the 
equations of gas dynamics use Rieraann solvers to resolve discontinuities at 
the interface between cells. A new method is proposed here for solving the 
Riemann problem based on a global existence proof for the solution to the 
Riemann problem. The method is found to be very reliable and computationally 
efficient. 
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1. INTRODUCTION 


The introduction of a finite-difference scheme based on Rieraann solvers by 
Godunov permitted the development of robust numerical methods for problems 
involving strong shocks. He proposed to consider all quantities in a 
computational cell as given by cell averages and to resolve the resulting 
discontinuities at the cell edges by the solution to a Riemann problem. The 
method has the important advantage that it automatically handles strong shocks 
and interactions and is able to predict cavitation should it occur. 

The use of Godunov's method has attracted particular interest in the field 
of aerodynamics. Here one expects to encounter only ideal gas (polytropic) 
equations of state. A Riemann solver for the equations of gas dynamics was 
first proposed by Godunov. It was subsequently modified and improved by 
Chorin [2], 

In this paper we take advantage of a global existence proof for the 
Riemann problem to propose an altogether different method for solving the 
Riemann problem. This method has the advantage that it is computationally 
very efficient and is guaranteed to converge to the solution of the Riemann 
problem. 


2. THE RIEMANN PROBLEM IN GAS DYNAMICS 

The gas dynamics equations in Eulerian coordinates are 


P t + (pu) x = 0 


(pu) fc + ( pu^ 


+ p) = 

r x 



e + ((e + p)u) = 0 

L X 


where p is the density of the gas, u is the velocity, pu is the momentum, 
e is the energy per unit volume, and p is the pressure. 

We assume that the gas is ideal and polytropic so that the equation of 
state is given by 

p = pRT 


where p is the pressure, p the density, T the temperature, and R the 
universal gas constant. Further, we assume that y > 1, where y is the 
ratio of specific heats. 

We seek a solution to the gas dynamics equations with initial data 



for 


x < 0 


and 



for 


x > 0. 


The solution to the Riemann problem is illustrated in Figure 1 on a 


x — t graph 
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Figure 1: Solution to the Riemann problem. 

The initial discontinuity is resolved by a system of waves consisting of a 
1-shock or rarefaction, a 2-contact discontinuity, and a 3-shock or rare- 
faction. If cavitation is excluded the pressure and velocity on either side 
of the contact discontinuity are equal and are given by p and u . 


3* THE GLOBAL EXISTENCE THEOREM 

A parametric representation of the solution to the Riemann problem and 
global existence proof are discussed in Smoller [3]. We outline the 
discussion to relate it to our numerical scheme. Let 


and 
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= y + 1 

Y “ 1 
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c L is the 


We denote the sound speed by c, where c = f— 1 ^ . Thus 
sound speed for the state The solution to the Rieraann problem is 

obtained by the following three one-parameter family of curves. 

1 - family, for Xj € 1R, 
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2 - family, for X 2 € 1 
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3 - family, for X3 £ IR 


P* 



Pp 

— = f 3 (x 3 ) 
p R 


f l (x 3 ) 
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U R - U 

^ h 3 ( x 3 ) 
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Y " 1 


1/9 X 3 

2(t) /2 (e - 1) 

Y ' ‘ (1 ♦ V 2 


x 3 ^ 0 


x 3 < 0. 


To obtain a solution of the Rieraann problem, we need to solve the above 
system of equations for the parameters x^ , X2, and X3. It should be noted 
that x^ < 0 corresponds to a 1-shock and x^ > 0 to a 1-rarefaction wave. 
A similar statement holds for the parameter X3. Using the above representa- 
tion the following theorem can be proved. 


Theorem: Let ( p L> u L ,p L^ an< * (p r ,Ur,Pr) be any two states (not 
necessarily close). Then there is a unique solution to the Riemann problem 
with these initial states if and only if 


U R " 


U T < — 
L Y 


- 1 


<C L + V- 


If this condition is violated, then a vacuum is present in the solution. 



4. THE NUMERICAL METHOD 


We first check whether 


u 


R 


- u T 


y - i 


(c T 


+ C R } - 


If the condition is violated, then cavitation occurs. 

To obtain the values x^ , X 2 , and X 3 explicitly, we define 



We first solve for x^, where x^ satisfies the equation: 
g(Xj) = hjCxj) + (^-) ^hjCxj + log B) - C = 0. 


It should be noted that the function g is a monotone function* Hence, the 
equation can be solved efficiently by a numerical routine. We first locate an 
interval [a,b] such that g(a) < 0 and g(b) > 0. With a and b as our 
initial guesses, a few iterations of the regula-falsi method are needed to 
obtain an accurate value for x^ . 

We then have 


x 3 = X 1 + log B. 
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Finally, we obtain X 2 from the equation: 

X 2 

fj (x 1 )e f 3^ x 3) = A * 

obtained , X 2 > and x^ we substitute them into the parametric 

representaion for the solution to the Rieraann problem and obtain the values 

* * * * 
u > P > P L > and p R . 


5. NUMERICAL RESULTS 

We found that over a range of y and on a variety of shock tube problems 
our Riemann solver converged to the exact solution within about five 
iterations of the regula-falsi method. 


6. CONCLUSIONS 

A numerical method for solving the Riemann problem in gas dynamics has 
been proposed. The method is based on an existence proof for the solution to 
the Riemann problem. Its advantages are that it can predict the occurrence of 
cavitation and it is computationally very efficient. 
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AP PEND IX 

Computer Program for the Riematm Solver 


subroutine riemann 
common/ input /ul ,ur, rhol, rhor, pi, pr 
common/ output /pstar , ustar 
common/parm/ gamma , beta , t au , a , b , c 
c input:ul,ur — velocities to the left and right of interface 

c rhol, rhor — densities 

c pl,pr — pressures 

c output: pstar — interface pressure 

c ustar — interface velocity 

beta=(gamma+l . )/ (gamma-1 . ) 
tau=( gamma- 1 . )/ (2.*gamma) 
cl=sqrt(gamma*pl/rhol) 
cr=sqrt(gamma*pr/rhor) 
c check for cavitaiton 

if ((ur-ul).gt.(2.*(cl+cr)/(gamma-l. )))go to 3 

a=rhor/rhol 

b=pr/pl 

c=(ur-ul)/cl 

x=0. 

c locate interval on which g(x) has a zero 
9 if ((g(x).gt .0. ).or.(g(x+l . ).lt .0. ) )then 
if (g(x).gt.0.)then 
x=x-l . 

else 

x=x+l . 
end if 

end if 

if ( (g(x) .gt .0 . ) • or. (g(x+l . ) .It .0. ) )go to 9 
c do regula-falsi iterations 
xmax=x+l . 
xmin=x 
n=l 

8 x=xmax-g(xmax)*(xmax-xmin)/ (g(xmax)-g(xmin) ) 
if (g(x).lt .0.) then 
xrain=x 

else 

xmax^x 
end if 

n=n+l 

if ( (abs(xmax-xmin) .gt . 0.0001). and. (n.le. 10)) go to 8 
xl=xraax 

c calculate ustar, pstar 
pstar=exp(-l .*xl )*pl 
ustar=ul+cl*hl (xl ) 
return 
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3 print *, 'cavitation occurs' 
return 
end 

function hl(x) 

common/ parm/ gamma , beta, tau , a , b, c 
if (x.gt.O. )then 

hi =( 2. / (gamma- 1 . ) )*(1 .-exp(-l . *tau*x) ) 

else 

hl=2.*sqrt(tau)*(l .-exp(-l .*x))/ (gamma-1 . ) 
hl=hl/sqrt(l .+beta*exp(-l .*x) ) 

end if 
return 
end 

function g(x) 
external hi 

common/parra/ gamma , beta , tau , a , b , c 
y=x+alog(b) 

g=hl(x)+sqrt(b/a)*hl (y)-c 

return 

end 
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